rm(list=ls(all.names=TRUE))
set.seed(221269)


#### Loading data 
load("anes2000.RData")


#### Loading libraries
library(lavaan)
library(semTools)



#### Setting ordered variables ####
orderedvars <- c('respect', 'manners', 'obedience', 'behaved')



####################
#### Configural ####
####################

Config1 <- 
  '
## Factor loadings
authorit =~ c(NA,NA)*respect
authorit =~ c(NA,NA)*manners
authorit =~ c(NA,NA)*obedience
authorit =~ c(NA,NA)*behaved


## Thresholds
respect | c(NA,NA)*t1
respect | c(NA,NA)*t2

manners | c(NA,NA)*t1
manners | c(NA,NA)*t2

obedience | c(NA,NA)*t1
obedience | c(NA,NA)*t2

behaved | c(NA,NA)*t1
behaved | c(NA,NA)*t2


## Observed variables (variate): Latent variate scale
respect     ~*~ c(1,1)*respect
manners     ~*~ c(1,1)*manners
obedience   ~*~ c(1,1)*obedience
behaved     ~*~ c(1,1)*behaved


#Latent variables: variances
authorit ~~ c(1,1)*authorit


## Observed variables (latent variate): Intercepts
respect    ~ c(0,0)*1
manners    ~ c(0,0)*1
obedience  ~ c(0,0)*1
behaved    ~ c(0,0)*1


#Latent variable: mean
authorit ~ c(0,0)*1


## Residual correlation, freely est
# respect ~~ c(NA,NA)*behaved
'

Config1fit <- lavaan::lavaan(
  data=anes2000, 
  model=Config1, 
  group="race", 
  group.label=c("White", "Black"), 
  ordered=orderedvars, 
  meanstructure=TRUE, parameterization="delta", 
  test="scaled.shifted", se="robust.sem", estimator="ULS",
  mimic="lavaan", model.type="cfa", start="simple", 
  # Weights
  sampling.weights="weight", 
  sampling.weights.normalization="group", 
  # Identification constraints
  int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
  std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
  auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
  auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
  # Categorical estimation
  zero.cell.warn=TRUE, 
  zero.add=c(0.5,0.5), 
  zero.keep.margins=TRUE, 
  # Computational / optimization options
  check.start=TRUE, check.post=TRUE,
  check.gradient=TRUE, check.vcov=TRUE, 
  check.lv.names=FALSE,
  do.fit=TRUE, optim.method="nlminb",
  implied=TRUE, verbose=TRUE)

summary(Config1fit, fit=T, standardized=T)
lavResiduals(Config1fit)
# fitMeasures(Config1fit)
modificationIndices(Config1fit)



#### Residual corr added ####

Config2 <- 
  '
## Factor loadings
authorit =~ c(NA,NA)*respect
authorit =~ c(NA,NA)*manners
authorit =~ c(NA,NA)*obedience
authorit =~ c(NA,NA)*behaved


## Thresholds
respect | c(NA,NA)*t1
respect | c(NA,NA)*t2

manners | c(NA,NA)*t1
manners | c(NA,NA)*t2

obedience | c(NA,NA)*t1
obedience | c(NA,NA)*t2

behaved | c(NA,NA)*t1
behaved | c(NA,NA)*t2


## Observed variables (variate): Latent variate scale
respect     ~*~ c(1,1)*respect
manners     ~*~ c(1,1)*manners
obedience   ~*~ c(1,1)*obedience
behaved     ~*~ c(1,1)*behaved


#Latent variables: variances
authorit ~~ c(1,1)*authorit


## Observed variables (latent variate): Intercepts
respect    ~ c(0,0)*1
manners    ~ c(0,0)*1
obedience  ~ c(0,0)*1
behaved    ~ c(0,0)*1


#Latent variable: mean
authorit ~ c(0,0)*1


## Residual correlation, freely est
respect ~~ c(NA,NA)*behaved
'

Config2fit <- lavaan::lavaan(
  data=anes2000, 
  model=Config2, 
  group="race", 
  group.label=c("White", "Black"), 
  ordered=orderedvars, 
  meanstructure=TRUE, parameterization="delta", 
  test="scaled.shifted", se="robust.sem", estimator="ULS",
  mimic="lavaan", model.type="cfa", start="simple", 
  # Weights
  sampling.weights="weight", 
  sampling.weights.normalization="group", 
  # Identification constraints
  int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
  std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
  auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
  auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
  # Categorical estimation
  zero.cell.warn=TRUE, 
  zero.add=c(0.5,0.5), 
  zero.keep.margins=TRUE, 
  # Computational / optimization options
  check.start=TRUE, check.post=TRUE,
  check.gradient=TRUE, check.vcov=TRUE, 
  check.lv.names=FALSE,
  do.fit=TRUE, optim.method="nlminb",
  implied=TRUE, verbose=TRUE)

summary(Config2fit, fit=T, standardized=T)
lavResiduals(Config2fit)
# fitMeasures(Config2fit)
modificationIndices(Config2fit)



##########################
#### Thresh + Loading ####
##########################

LoadThresh1 <- 
  '
## Factor loadings
authorit =~ c(LoRespWh,LoRespBl)*respect 
authorit =~ c(LoManrWh,LoManrBl)*manners 
authorit =~ c(LoObedWh,LoObedBl)*obedience 
authorit =~ c(LoBehaWh,LoBehaBl)*behaved
LoRespWh == LoRespBl
LoManrWh == LoManrBl
LoObedWh == LoObedBl
LoBehaWh == LoBehaBl


## Thresholds
respect | c(ThRespWh1,ThRespBl1)*t1
respect | c(ThRespWh2,ThRespBl2)*t2
ThRespWh1 == ThRespBl1
ThRespWh2 == ThRespBl2

manners | c(ThManrWh1,ThManrBl1)*t1
manners | c(ThManrWh2,ThManrBl2)*t2
ThManrWh1 == ThManrBl1
ThManrWh2 == ThManrBl2

obedience | c(ThObedWh1,ThObedBl1)*t1
obedience | c(ThObedWh2,ThObedBl2)*t2
ThObedWh1 == ThObedBl1
ThObedWh2 == ThObedBl2

behaved | c(ThBehaWh1,ThBehaBl1)*t1
behaved | c(ThBehaWh2,ThBehaBl2)*t2
ThBehaWh1 == ThBehaBl1
ThBehaWh2 == ThBehaBl2


## Observed variables (variate): Latent variate scale
respect     ~*~ c(1,NA)*respect
manners     ~*~ c(1,NA)*manners
obedience   ~*~ c(1,NA)*obedience
behaved     ~*~ c(1,NA)*behaved


#Latent variables: variances
authorit ~~ c(1,NA)*authorit


## Observed variables (latent variate): Intercepts
respect    ~ c(0,NA)*1
manners    ~ c(0,NA)*1
obedience  ~ c(0,NA)*1
behaved    ~ c(0,NA)*1


#Latent variable: mean
authorit ~ c(0,0)*1


## Residual correlation, freely est
respect ~~ c(ResWh,ResBl)*behaved
# ResWh == ResBl
'

LoadThresh1fit <- lavaan::lavaan(
  data=anes2000, 
  model=LoadThresh1, 
  group="race", 
  group.label=c("White", "Black"), 
  ordered=orderedvars, 
  meanstructure=TRUE, parameterization="delta", 
  test="scaled.shifted", se="robust.sem", estimator="ULS",
  mimic="lavaan", model.type="cfa", start="simple", 
  # Weights
  sampling.weights="weight",
  sampling.weights.normalization="group",
  # Identification constraints
  int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
  std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
  auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
  auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
  # Categorical estimation
  zero.cell.warn=TRUE, 
  zero.add=c(0.5,0.5), 
  zero.keep.margins=TRUE, 
  # Computational / optimization options
  check.start=TRUE, check.post=TRUE,
  check.gradient=TRUE, check.vcov=TRUE, 
  check.lv.names=FALSE,
  do.fit=TRUE, optim.method="nlminb",
  implied=TRUE, verbose=TRUE)

summary(LoadThresh1fit, fit=T, standardized=T)
lavResiduals(LoadThresh1fit)
# fitMeasures(LoadThresh1fit)
lavTestScore(LoadThresh1fit, univariate=TRUE, epc=TRUE)



###################
#### Intercept ####
###################

Intercept1 <- 
'
## Factor loadings
authorit =~ c(LoRespWh,LoRespBl)*respect 
authorit =~ c(LoManrWh,LoManrBl)*manners 
authorit =~ c(LoObedWh,LoObedBl)*obedience 
authorit =~ c(LoBehaWh,LoBehaBl)*behaved
LoRespWh == LoRespBl
LoManrWh == LoManrBl
LoObedWh == LoObedBl
LoBehaWh == LoBehaBl


## Thresholds
respect | c(ThRespWh1,ThRespBl1)*t1
respect | c(ThRespWh2,ThRespBl2)*t2
ThRespWh1 == ThRespBl1
ThRespWh2 == ThRespBl2

manners | c(ThManrWh1,ThManrBl1)*t1
manners | c(ThManrWh2,ThManrBl2)*t2
ThManrWh1 == ThManrBl1
ThManrWh2 == ThManrBl2

obedience | c(ThObedWh1,ThObedBl1)*t1
obedience | c(ThObedWh2,ThObedBl2)*t2
ThObedWh1 == ThObedBl1
ThObedWh2 == ThObedBl2

behaved | c(ThBehaWh1,ThBehaBl1)*t1
behaved | c(ThBehaWh2,ThBehaBl2)*t2
ThBehaWh1 == ThBehaBl1
ThBehaWh2 == ThBehaBl2


## Observed variables (variate): Latent variate scale
respect     ~*~ c(1,NA)*respect
manners     ~*~ c(1,NA)*manners
obedience   ~*~ c(1,NA)*obedience
behaved     ~*~ c(1,NA)*behaved


#Latent variables: variances
authorit ~~ c(1,NA)*authorit


## Observed variables (latent variate): Intercepts
respect    ~ c(IntRespWh,IntRespBl)*1
manners    ~ c(IntManrWh,IntManrBl)*1
obedience  ~ c(IntObedWh,IntObedBl)*1
behaved    ~ c(IntBehaWh,IntBehaBl)*1
IntRespWh == 0 
IntRespWh == IntRespBl
IntManrWh == 0
IntManrWh == IntManrBl
IntObedWh == 0
IntObedWh == IntObedBl
IntBehaWh == 0
IntBehaWh == IntBehaBl


#Latent variable: mean
authorit ~ c(0,NA)*1


## Residual correlation, freely est
respect ~~ c(ResWh,ResBl)*behaved
# ResWh == ResBl
'

Intercept1fit <- lavaan::lavaan(
  data=anes2000, 
  model=Intercept1, 
  group="race", 
  group.label=c("White", "Black"), 
  ordered=orderedvars, 
  meanstructure=TRUE, parameterization="delta", 
  test="scaled.shifted", se="robust.sem", estimator="ULS",
  mimic="lavaan", model.type="cfa", start="simple", 
  # Weights
  sampling.weights="weight", 
  sampling.weights.normalization="group", 
  # Identification constraints
  int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
  std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
  auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
  auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
  # Categorical estimation
  zero.cell.warn=TRUE, 
  zero.add=c(0.5,0.5), 
  zero.keep.margins=TRUE, 
  # Computational / optimization options
  check.start=TRUE, check.post=TRUE,
  check.gradient=TRUE, check.vcov=TRUE, 
  check.lv.names=FALSE,
  do.fit=TRUE, optim.method="nlminb",
  implied=TRUE, verbose=TRUE)

summary(Intercept1fit, fit=T, standardized=T)
lavResiduals(Intercept1fit)
# fitMeasures(Intercept1fit)
(lav_test_int <- lavTestScore(Intercept1fit, univariate=TRUE, epc=TRUE))

restrict <- which(lav_test_int$uni$rhs=="IntRespBl")
epc_IntRespBl <- lavTestScore(Intercept1fit, release=restrict, 
                              univariate=TRUE, epc=TRUE)
(epc_interest <- subset(epc_IntRespBl$epc, op=="~1"))


# Release intercept RESPECT
Intercept2 <- 
  '
## Factor loadings
authorit =~ c(LoRespWh,LoRespBl)*respect 
authorit =~ c(LoManrWh,LoManrBl)*manners 
authorit =~ c(LoObedWh,LoObedBl)*obedience 
authorit =~ c(LoBehaWh,LoBehaBl)*behaved
LoRespWh == LoRespBl
LoManrWh == LoManrBl
LoObedWh == LoObedBl
LoBehaWh == LoBehaBl


## Thresholds
respect | c(ThRespWh1,ThRespBl1)*t1
respect | c(ThRespWh2,ThRespBl2)*t2
ThRespWh1 == ThRespBl1
ThRespWh2 == ThRespBl2

manners | c(ThManrWh1,ThManrBl1)*t1
manners | c(ThManrWh2,ThManrBl2)*t2
ThManrWh1 == ThManrBl1
ThManrWh2 == ThManrBl2

obedience | c(ThObedWh1,ThObedBl1)*t1
obedience | c(ThObedWh2,ThObedBl2)*t2
ThObedWh1 == ThObedBl1
ThObedWh2 == ThObedBl2

behaved | c(ThBehaWh1,ThBehaBl1)*t1
behaved | c(ThBehaWh2,ThBehaBl2)*t2
ThBehaWh1 == ThBehaBl1
ThBehaWh2 == ThBehaBl2


## Observed variables (variate): Latent variate scale
respect     ~*~ c(1,NA)*respect
manners     ~*~ c(1,NA)*manners
obedience   ~*~ c(1,NA)*obedience
behaved     ~*~ c(1,NA)*behaved


#Latent variables: variances
authorit ~~ c(1,NA)*authorit


## Observed variables (latent variate): Intercepts
respect    ~ c(IntRespWh,IntRespBl)*1
manners    ~ c(IntManrWh,IntManrBl)*1
obedience  ~ c(IntObedWh,IntObedBl)*1
behaved    ~ c(IntBehaWh,IntBehaBl)*1
IntRespWh == 0 
# IntRespWh == IntRespBl
IntManrWh == 0
IntManrWh == IntManrBl
IntObedWh == 0
IntObedWh == IntObedBl
IntBehaWh == 0
IntBehaWh == IntBehaBl


#Latent variable: mean
authorit ~ c(0,NA)*1


## Residual correlation, freely est
respect ~~ c(ResWh,ResBl)*behaved
# ResWh == ResBl
'

Intercept2fit <- lavaan::lavaan(
  data=anes2000, 
  model=Intercept2, 
  group="race", 
  group.label=c("White", "Black"), 
  ordered=orderedvars, 
  meanstructure=TRUE, parameterization="delta", 
  test="scaled.shifted", se="robust.sem", estimator="ULS",
  mimic="lavaan", model.type="cfa", start="simple", 
  # Weights
  sampling.weights="weight", 
  sampling.weights.normalization="group", 
  # Identification constraints
  int.ov.free=TRUE, int.lv.free=TRUE, auto.fix.first=FALSE,
  std.lv=FALSE, orthogonal=FALSE, auto.fix.single=FALSE,
  auto.th=TRUE, auto.delta=TRUE, auto.var=TRUE, 
  auto.cov.lv.x=TRUE, auto.cov.y=TRUE,
  # Categorical estimation
  zero.cell.warn=TRUE, 
  zero.add=c(0.5,0.5), 
  zero.keep.margins=TRUE, 
  # Computational / optimization options
  check.start=TRUE, check.post=TRUE,
  check.gradient=TRUE, check.vcov=TRUE, 
  check.lv.names=FALSE,
  do.fit=TRUE, optim.method="nlminb",
  implied=TRUE, verbose=TRUE)

summary(Intercept2fit, fit=T, standardized=T)
lavResiduals(Intercept2fit)
# fitMeasures(Intercept2fit)
lavTestScore(Intercept2fit, univariate=TRUE, epc=TRUE)


# Comparing nested models: Full v. partial intercept
lavTestLRT(Intercept1fit, Intercept2fit)
